Power-law rheology and the dynamical heterogeneity in a sheared granular material 
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Rheology of a granular material at the jamming density is investigated using molecular dynamics 
simulation. It is found that shear stress exhibits power-law dependence on shear rate with a nontriv- 
ial exponent. Due to the criticality of the jamming transition point, finite-size effect is observed in 
smaller systems at lower shear rates. Finite-size scaling indicates the correlation length algebraically 
diverges in the zero shear rate limit. It is also found that the dynamical susceptibility monotoni- 
cally decreases with time so that the dynamical heterogeneity is detected by a two-point correlation 
function. Several exponents that describe rheology, the correlation length, and the amplitude of the 
dynamical susceptibility are estimated. 



Consider an assembly of repulsive particles at zero tem- 
perature, the density of which is so low that there is no 
contact between particles. By increasing density, the par- 
ticles are inevitably in contact above a certain density 
and the system obtains rigidity without crystallization 
P, Q • This phenomenon and the threshold density are 
now referred to as the jamming transition and the jam- 
ming density, respectively. It is believed that the jam- 
ming transition point, which is also referred to as point 
J, is a critical point at zero temperature due to several 
evidences: finite-size scaling of the transition point 0, 
and the growing correlation length and relaxation time 

0, 0, SH, 0, B However, in contrast to conven- 
tional critical phenomena, these length and time scales 
are apparent only through dynamic quantities such as the 
displacement of particles in a certain lag time. Quite in- 
terestingly, this dynamic critical fluctuation is originally 
observed in supercooled liquids 13, 111, and is referred to 
as the dynamical heterogeneity. Recalling that the jam- 
ming transition is the solidification without crystalliza- 
tion, the jamming and glass transitions might be closely 
related. 

However, the nature of point J is still opaque be- 
cause the critical exponents and the mechanism that 
yields the dynamical heterogeneity are not clear. As re- 
gards the critical exponents, it is generally accepted that 
C ~ \4>J — 4'\~'^ : where (j) and denote the density of 
a system and the jamming density. The critical expo- 
nent 1/ is estimated as 0.6 ± 0.1 d, H, [1, S [13] ■ However, 
the exponents for other relevant variables are not known. 
For example, because point J is defined at zero temper- 
ature and zero shear rate limit P, [3], temperature T 
and shear rate 7 are also relevant variables so that the 
correlation length must diverge as T ^ and 7^0. 
However, we still do not know the nature of divergence; 

1. e., power law, Vogel-Fulcher law, etc. Unfortunately, 
the concept of temperature in granular media is not gen- 
erally established at this point so that we focus only on 
shear rate. 

In this paper, by molecular dynamics simulation, we 
investigate rheology of a granular material at the jam- 
ming density focusing on the particle dynamics and the 



correlation length. Several exponents that describe rhe- 
ology, the correlation length, and the amplitude of the 
dynamical susceptibility are estimated. 

Our simulation system consists of monodisperse 
spheres, the diameter of which is denoted by d. The 
repulsive force between particles is given as k{d — |r|) if 
d > |r| and otherwise zero, where r denote the relative 
vector of the position of two particles. The inelasticity is 
modeled by (Cn • r)n, where n — r/|r|. The total inter- 
action is the sum of the above two terms. Throughout 
this study, we adopt the units in which d = 1, k = 1, 
and m = 1, where m is the mass of a particle. We set 
C — 2, which corresponds to the damping limit; i.e., the 
coefficient of restitution is zero. 

The dimensions of the system is L x L x L, where L 
is varied from 3.45 to 40.32. Here we fix the density of 
the system to be 0.639 in the volume density, which is 
considered as the jamming density for three dimensional 
systems with the precision of ±0.001. The largest system 
{L — 40.32) contains 80, 000 particles, while the smallest 
system consists of only 50 particles. We set the flow di- 
rection along the y axis and the velocity gradient along 
the z axis, and adopt Lees-Edwards boundary conditions 
[1^ . The equations of motion are given by SLLOD equa- 
tions; 



= 

Pi = 



\- iqz.t^y, 

m 

Fint - IPz.i'n.y, 



(1) 

(2) 



where Fint denotes the repulsive interaction force be- 
tween particles and liy is the unit vector along the y axis. 
In the present configuration, the control parameters are 
the shear rate 7 and the system size L so that the shear 
stress S and the pressure P are formally expressed as 
S = 5(7, L-i) and P = F(7, L'^). 

First we discuss power-law rheology shown in FIG. [1] 
In the higher shear rate region, the pressure and the shear 
stress exhibit power-law dependence on shear rate, each 
of which has different exponents: S oc 'j^^^'^ and P oc 
7!/-^^, where l/Ss = 0.64 ±0.02 and l/5p = 0.49 ±0.01. 
Note that these exponents differ from 2.0 (i.e., Bagnold's 
scaling ^sj), which is valid only in imjammed regime 
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FIG. 1: System-size dependence of the rheological properties, 
(a) Shear-rate dependence of shear stress. The dashed line 



denotes 5* oc 7 



(b) Shear-rate dependence of pressure. 



The dashed line denotes P oc 7 
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FIG. 2: Collapse of the rheological data shown in FIG. [T] by 
finite-size scaling: (a) Eq. (|3]), and (b) Eq. Q. 



where the stiffness of the particles does not affect the 
mean free time. Note that this power-law rheology with 
nontrivial exponent is one of the evidences of the criti- 
cality of point J. We also remark that these exponents 
do not contradict those previously conjectured [li 



These power-law behaviors, however, do not hold in 
the lower shear rate region and there exists crossover. 
Moreover, the shear rate at which the crossover occurs 
depends on the system size; it shifts lower as the system 
becomes larger. Indeed, this crossover is described by 
finite-size scaling laws, which are of the same form as 
those in conventional critical phenomena [l8| . 



(3) 
(4) 



where fs and fp are scaling functions. FIG. [2] shows the 
collapse of the rheological data using Eqs. ([3|) and ([4]), 
where the exponents are estimated as 1/ys — 2.5 ± 0.2, 
l/yp = 2.0 ± 0.2, and l/y^ = 4.0 ± 0.5. 

We remark three important properties that are de- 
duced from Eqs. ^ and First, the exponents for 



power-law rheology, 5$ and 5p, are expressed by y^, ys, 
and yp. Because the rheology does not depend on the 
system size in the thermodynamic limit L — > 00, the scal- 
ing functions must behave as fs{z) cx z^i/y^ and fp{z) cx 
zVi/yp ^ respectively. They lead to £'(7,00) oc ^y-i/y^ and 
P(7, 00) cx ^y-i/yp. We thus obtain equalities for the 
exponents, l/5s — y-y/ys and 1/Sp = y-y/yp, which are 
confirmed by inserting the above values. Second, because 
finite-size crossover occurs at L/^^y'i 1, we can expect 
that there exist the growing correlation length that de- 
pends on the shear rate as 



(5) 



where y^ = 0.25 ± 0.03. We show that, in the following 
paragraph, this growing correlation length is explicitly 
observed in a spatial correlation function. Third, insert- 
ing L = ^ into Eqs. jH]) and (gl), we obtain 5' ~ ^"i/ws 
and P ^ ^-^/yp ^ which represents the scaling dimensions 
of the shear stress and the pressure, respectively. 

From Eq. (I5|), it is expected that the dynamical het- 
erogeneity is enhanced in the 7 — > limit. In order to 
detect the dynamical heterogeneity, we define the four- 
point correlation function h{r,T) = (/i((r,T))t, where 



ht(r,T) = ^^6m^{T;t)6mj{T;t)S{r ~ Xi{t) +Xj{t)). 

(6) 



Here Smiij^t) = mi(T;t) — {mi(T]t))i, where 



mi{T;t) = 
5xi{T;t) = 



exp 



SxjjT; t) 



m 



(7) 
(8) 



Using /i(r, r), we introduce the dynamical susceptibility 
X4(t) = J drh{r,T). As is shown in FIG. [3l the ampli- 
tude of X4(''') can be scaled by ^(7), which is proportional 
to 7"^^. The exponent is estimated as c = 0.65 ± 0.05. 
More importantly, this dynamical susceptibility seems to 
be a monotonically decreasing function with the lag time 
r. This means that the dynamical heterogeneity can be 
detected in the r — > limit; i.e., two-point correlation 
of the instantaneous velocity field. This makes an appar- 
ent contrast to supercooled liquids and vibrated granular 
matters, in which the dynamical susceptibility reaches 
maximum after a considerable lag time. This is because 
the particles are continuously rearranged due to exter- 
nally applied shear so that the agitation without struc- 
tural relaxation is absent in the present system. In addi- 
tion, the thermal motion of particles, which predominates 
the structural relaxation, is absent in granular matters. 
(In particular, the coefficient of restitution is zero in the 
present system.) Due to these points, the instantaneous 
velocity leads to structural relaxation. 

In the rest of this paper, we thus investigate the nature 
of the instantaneous velocity field. Figure 2] shows the 
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FIG. 3: The temporal behavior of the dynamical susceptibil- 
ity Xii^)- The lag time r is multiplied by the shear rate 7 
so that the horizontal axis represents strain. Note that the 
dynamical susceptibility X^i''') is normalized by ^(7). (Inset) 
The amplitude ^(7) is proportional to 7""'^^, shown as the 
solid line. 




FIG. 4: The correlation function defined by Eqs. ((SJ and ((9}. 
The legends are the same as those in FIG. |3l Note that the 
length scale is normalized by ^. (Inset) The correlation length 
^ is proportional to 7"'''^^, shown as the solid line. 



two-point correlation function defined in the r — > limit 
ofEq. ®. 

Hir)^ ^ f dnh{r,0), (9) 
g(r)h{0,0) J 

where g(r) denotes the radial distribution function and ft 
denotes the solid angle. This correlation function clearly 
shows that the correlation length ^ is proportional to 
j'y-' , where = 0.23 ± 0.02. Note that this is consis- 
tent with Eq. ([5|), which results from finite-size scaling. 
Therefore we can conclude that the finite-size effect in 
rheology occurs as the correlation length is comparable 
to the system size. 

Note that H{r) defined by Eqs. ([6|) and ^ character- 
izes the spatial heterogeneity of mobility, but does not 
have information on the specific motion of the particles 
of high mobility. In the kinetic theoretical point of view, 




FIG. 5: The correlation function defined by Eq. (I10|l . which 
represents the extent of stringlike cooperative motion. Note 
that the length scale is normalized by the characteristic length 
A. The legends are the same as those in FIG. |31 (Inset) The 
characteristic length A(7) is proportional to 7""'^^, shown as 
the solid line. 

it is important to know how particles rearrange cooper- 
atively in a marginally-jammed system. To this end, we 
define a spatial correlation function for the instantaneous 
velocity vectors; 

^(^'*) ^ ^5]p.W •p,W<5(r-x,(t)+x,(t)|). (10) 

This characterizes the spatial correlation of co-moving 
particles, such as ring collision in a dilute gas and 
the stringlike cooperative motion (SCM) in glassy sys- 
tems 0, [1^ • In addition, note that particles with larger 
velocity predominate G(r, t), because the magnitude of 
velocity is taken into account in Eq. PH)) . Here we 
average the correlation function with respect to time 
G{r) = {G{r,t))t, and the solid-angle dependence is in- 
tegrated out as is done in Eq. ([9]). After this procedure, 
the correlation function G{r) becomes an exponentially 
decaying function so that it can be scaled by a length 
scale A, which should be interpreted as the characteristic 
length of SCM. The scaled correlation function is shown 
in FIG. [5l The characteristic length of SCM increases 
as the shear rate decreases as is shown in the inset of 
FIG. [H The exponent is estimated as 0.15 if we assume 
power law; i.e., A ^-o i5^ This is to be compared with 
Eq. dni), from which we can conclude that SCM is not 
directly related to the dynamical heterogeneity. This sit- 
uation is similar in two-dimensional air-fluidized beads 
, in which several characteristic lengths for SCM are 
found. The kinetics of the dynamical heterogeneity is 
thus still open. 

The discussions so far constitute the main results of 
this study. In the following, we remark two points that 
are peripherally related to the main results. The first 
point is the comparison of the present result with Ref. 
9], in which a scaling law for the correlation length is 
found for a two-dimensional overdamp system. This im- 
plies that the correlation length diverges as ^ oc S~^^^. 
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TABLE I: The exponents estimated in the present study. 



Observable 



Critical 
exponents 



Shear stress Pressure Length Susceptibility 



l/ys l/Ss l/yp l/Sp 
2.5(2) 0.64(2) 2.0(2) 0.49(1) 0.23(2) 



0.65(5) 



Because they obtain S oc 7°-*^ at the critical density, 
this leads to ^ ~ 7"°-^^. The latter is consistent to the 
present result, while the former is different. The reason 
for this discrepancy is not due to the dimensionality but 
the artificial overdamp (massless) dynamics in Ref. 
because we obtain S oc 7*''^^ in a two-dimensional gran- 
ular system and also S oc 7^'-^^ in a three-dimensional 
overdamp system, the details of which are to be presented 
elsewhere. 

Then we wish to remark the experimental validation 
of Eq. ([5]) with ~ 0.25. Although this exponent is so 
small that the correlation length shown in FIG. H] does 
not significantly increase (by approximately one decade) , 
it is observed to increase by two order of magnitudes in 
an experiment [2l| |. 

In summary, power-law rheology is found in a granu- 
lar material at the jamming density. Finite-size scalings 
of power-law rheology indicate the diverging correlation 
length in the zero shear rate limit. The dynamical suscep- 
tibility becomes maximum for instantaneous velocity so 
that the correlation length of the dynamical heterogene- 
ity is defined by a two-point correlation function. This 
correlation length directly supports the result of finite- 
size scaling. Various critical exponents are estimated, as 
summarized in Table. [H 
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